
function [r] = R_Solve(q)
global  h   dx x  ;

dx2 = 1/(dx*dx);

%---- Construct the spatial Discretization Matrix of R------%

n = -2/15*h*h*h*dx2;
E=zeros(x,x);
E(1,x)= n;
E(1,2)= n ;
E(1,1)= -2*n+h/3;

for k=2:x-1
    E(k,k-1)=n;
    E(k,k+1) =n;
    E (k,k) = -2*n+h/3;
end

E(x,x-1)= n;
E(x,1)= n;
E(x,x)= -2*n+h/3;

b=h*h/3*dx2;

C=zeros(x,x);
C(1,x)= b;
C(1,2)= b;
C(1,1)= -2*b;
for k=2:x-1
    C(k,k-1)=b;
    C(k,k+1) =b;
    C (k,k) = -2*b;
end
C(x,x-1)= b;
C(x,1)= b;
C(x,x)= -2*b;

% Solve Er=-Dq
r = -E\(C*q);
